
restore,'./local/BGPS_V1_emaf.sav',/ver
restore,'./emaf/irdc_full_data.sav',/ver
restore,'./ancillary/irdc_s06_data.sav',/ver


help,pf_irdc,/str
help,s_irdc,/str

window,0,xsize=1000
multiplot,[2,1],xgap=0.035
plothist,pf_irdc.maj*2.,bin=10,charsize=1.0,xtitle='Major axis [arcsec]',$
         color='opposite'
plothist,s_irdc.majaxis*60.,bin=10,/over,color='red'

multiplot

plothist,pf_irdc.min*2.,bin=10,charsize=1.0,xtitle='Minor axis [arcsec]',$
         color='opposite'
plothist,s_irdc.minaxis*60.,bin=10,/over,color='red'

multiplot,/reset
multiplot,/default
multiplot,/reset

n = n_elements(emaf)
n_pf  = intarr(n)
n_s   = intarr(n)
n_any = intarr(n)

info = {l:0.,b:0.,c:0.,smaj:0.,smin:0.,dis:0.}

emaf_irdc = replicate({pf:info,s:info},n)

FOR i=0,n-1 DO BEGIN
   
   IF (i+1) MOD 100 EQ 0 THEN message,'Working...',/inf
   
   ;; Get center-to-center separation
   gcirc,2,emaf[i].l,emaf[i].b,pf_irdc.l,pf_irdc.b,pf_dis
   gcirc,2,emaf[i].l,emaf[i].b,s_irdc.glon,s_irdc.glat,s_dis
   
   ;; Are distances within circle of semi-major axis?
   is_pf = pf_dis LE pf_irdc.maj
   is_s  = s_dis  LE s_irdc.majaxis
   
   n_pf[i]  = fix(total(is_pf))
   n_s[i]   = fix(total(is_s))
   n_any[i] = n_pf[i] + n_s[i]
   
   ;; Load up the EMAF_IRDC structure
   IF n_pf[i] THEN BEGIN
      mindis = min(pf_dis,j)
      IF mindis LE 30. THEN BEGIN
         emaf_irdc[i].pf.l    = pf_irdc[j].l
         emaf_irdc[i].pf.b    = pf_irdc[j].b
         emaf_irdc[i].pf.smaj = pf_irdc[j].maj
         emaf_irdc[i].pf.smin = pf_irdc[j].min
         emaf_irdc[i].pf.c    = 1. - (pf_irdc[j].i_min/pf_irdc[j].i_mir)
         emaf_irdc[i].pf.dis  = mindis
      ENDIF
   ENDIF
   IF n_s[i] THEN BEGIN
      mindis = min(s_dis,k)
      IF mindis LE 30. THEN BEGIN
         emaf_irdc[i].s.l    = s_irdc[k].glon
         emaf_irdc[i].s.b    = s_irdc[k].glat
         emaf_irdc[i].s.smaj = s_irdc[k].majaxis
         emaf_irdc[i].s.smin = s_irdc[k].minaxis
         emaf_irdc[i].s.c    = s_irdc[k].peak
         emaf_irdc[i].s.dis  = mindis
      ENDIF
   ENDIF   
   
ENDFOR

save,emaf_irdc,filename='./emaf_paper/save_files/emaf_irdc.sav',/ver

myps,'./emaf_paper/plots/irdc_histograms.eps'

plothist,n_pf,color='opposite',xr=[-1,5],thick=3
plothist,n_s,/over,color='red',thick=3
plothist,n_any,/over,color='cyan',thick=3

print,'Ns: ',m4_stat(n_pf),m4_stat(n_s),m4_stat(n_any)

iz = where(n_any EQ 0, nz)
al_legend,box=0,/top,/right,$
          ['N!dzero!n = '+string(nz,format="(I0)"),$
           'Fraction = '+string(float(nz)/float(n),format="(F0.2)"),$
           'PF09','S06','Any'],linestyle=0,linsize=0.5,$
          color=[replicate('white',2),'opposite','red','cyan'],thick=3

          
myps,/done

myps,'./emaf_paper/plots/emaf_irdc_comps.eps',xsize=15
multiplot,[3,1]

cgPlot,charsize=1.0,n_pf,emaf.c,psym=16,symsize=0.5,xtit='N!dPF09!n',ytit='C',$
       xr=[-1,5]

linfit_plot,n_pf,emaf.c,color='cyan',/top,/left,/legend,linestyle=2

multiplot

cgPlot,charsize=1.0,n_s,emaf.c,psym=16,symsize=0.5,xtit='N!dS06!n',xr=[-1,5]
linfit_plot,n_s,emaf.c,color='cyan',/top,/left,/legend,linestyle=2

multiplot

cgPlot,charsize=1.0,n_any,emaf.c,psym=16,symsize=0.5,xtit='N!dIRDC!n',xr=[-1,5]
linfit_plot,n_any,emaf.c,color='cyan',/top,/left,/legend,linestyle=2

myps,/done,/mp
END
